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Abstract 

We present a sequential Monte Carlo sampler variant of the partial rejection 
control algorithm, and show that this variant can be considered as a sequen- 
tial Monte Carlo sampler with a modified mutation kernel. We prove that the 
new sampler can reduce the variance of the incremental importance weights 
when compared with standard sequential Monte Carlo samplers. We provide 
a study of theoretical properties of the new algorithm, and make connections 
with some existing algorithms. Finally, the sampler is adapted for applica- 
tion under the challenging "likelihood free," approximate Bayesian computa- 
tion modelling framework, where we demonstrate superior performance over 
existing likelihood-free samplers. 
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1 Introduction 



Sequential Monte Carlo (SMC) methods have emerged out of the fields of engineer- 
ing, probability and statistics in recent years. Variants of the methods sometimes 
appear under the names of particle filtering or interacting particle systems (e.g. 
Arulampalam et al. 2002] lAndrieu et al. 20031 [Del Moral 20041 IDoucet et al. 2001[). 



and their theoretical properties have been extensively studied (ICrisan and Doucet~2002l 
IDel Moral 2004] IKunschl2005l) . 

The standard SMC algorithm involves finding a numerical solution to a set of 
filtering recursions, such as filtering problems arising from non-linear / non-Guassian 
state space models. Under this framework, the SMC algorithm samples from a (often 
naturally occurring) sequence of distributions ir t , indexed by t = 1,...,T. Each 



distribution is defined on the support E l = ExEx...xE. Del Moral et al. (2006) 



sec 



also lPeters 20051) generalize the SMC algorithm to the case where the distributions irt 
are all defined on the same support E. This generalization, termed the SMC sampler, 
adapts the SMC algorithm to the more popular setting in which the state space E 
remains static. 

In short, the SMC sampler generates weighted samples (termed particles) from 
a sequence of distributions irt, for t = 1, . . . ,T, where ir? may be of particular in- 
terest. We refer to ir? as the target distribution. Procedurally, particles obtained 
from an arbitrary initial distribution tt\, with a set of corresponding initial weights, 
are sequentially propagated through each distribution ir t in the sequence via three 
processes, involving mutation (or move), correction (or importance weighting) and 
selection (or resampling). The final weighted particles at distribution ttt are con- 
sidered weighted samples from the target distribution it. The mechanism is similar 
to sequential importance sampling (resampling) (ILiu 2001UDoucet et al. 20011) . with 
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one of the crucial differences being the framework under which the particles are al- 
lowed to move, resulting in differences in the calculation of the weights of the particles. 

One of the major difficulties with SMC-type algorithms is particle depletion, in 
which the weights of the majority of the particles gradually decrease to zero, while 
a few particle weights dominate the population. This severely increases the vari- 
ability of Monte Carlo estimates of expectations under ir. In this article, we de- 
velop an algorithm which incorporates the partial rejection control (PRC) strategy of 



Liu (2001) into the SMC sampler framework. A particular motivation for this stems 



from the recent developments in "likelihood-free" (or approximate Bayesian) com- 
putation (IBeaumont et al. 2002) Marjoram et al. 2003 ISisson et al. 2007)) . where an 
extremely high proportion of mutated particles are expected to have very small, or 
exactly zero, posterior weights. 

In this article, we develop the SMC samplers PRC algorithm, in which the partial 
rejection control mechanism is built directly into the mutation kernel of the SMC 
sampler. In this manner, a particle mutation may be rejected if the resulting impor- 
tance weight is below a certain threshold. We begin with a brief introduction to the 
standard sequential Monte Carlo sampler in Section [2} and then present the SMC 
sampler PRC algorithm. We also discuss implementational issues arising from the 
inclusion of the PRC stage, including estimation for the resultant kernel normalizing 
constant. Section [3] provides some theoretical results that justify the addition of PRC 
in terms of improvements in the variance of the incremental importance weights. We 
also discuss a central limit theorem and derive a recursive expression for the asymp- 
totic variance of our algorithm. In addition, we make a novel connection between the 
SMC sampler PRC algorithm and the AliveSMC algorithm from the rare-event lit- 



erature developed in Le Gland and Oudjane (2004) In Section 0] we adapt the SMC 
sampler PRC algorithm for application in the likelihood-free modelling framework, 
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and demonstrate the computational gains achieved over existing likelihood-free al- 
gorithms via a simulated example. Finally, we present a stochastic claims reserving 
analysis using the developed methods in Section and conclude with a discussion. 



2 Sequential Monte Carlo and partial rejection 
2.1 Sequential Monte Carlo sampler 



Del Moral et al. (2006) introduced a modification of the sequential Monte Carlo al- 
gorithm, termed the sequential Monte Carlo sampler. Consider a sequence of distri- 
butions ir t (x),t = 1, . . . ,T, with x 6 E, where the final distribution ttt is the distri- 
bution of interest. By introducing a sequence of backward kernels Lk, a new distribu- 
te 

tion 7r t (xi, . . . , x t ) = nt(%t) EI Lk (xk+i,Xk) may be defined for the path of a particle 

k=l 

(xi, . . . , x t ) £ E 1 through the sequence 7Ti, . . . , ~k%. The only restriction on the back- 
ward kernels is that the correct marginal distributions J TT t (xi, . . . , Xt)dxi, . . . , dxt-i = 
7Tt(xt) are available. 

Within this framework, one may then work with the sequence of distributions, 
TT t , under the standard SMC algorithm. In summary, the SMC sampler algorithm 
involves three stages: mutation, whereby the particles are moved from x t ~\ to x t 
via a mutation kernel M t (xt-i,Xt) as described below (IDoucet and Johansen 2009 



IDel Moral et al. 20061) ; correction, where the particles are reweighted with respect to 



ir t via the incremental importance weight ((T]); and selection, where according to some 
measure of particle diversity, commonly the effective sample size (ESS, IKunsch 2005( 
Kitigawa 1996| IDoucet et al. 20011 ILiu and Chen 19981) . the weighted particles may 



be resampled in order to reduce the variability of the importance weights. 

In more detail, suppose that at time t — 1, the distribution 7r t _i can be approxi- 
mated empirically by tt^Li using N weighted particles. These particles are first prop- 



agated to the next distribution TT t using a mutation kernel M t (x t -i,x t ), and then 
assigned new weights W t = W t -iw t (xi, . . . x t ), where Wt-i is the weight of a particle 
at time t — 1 and Wt is the incremental weight given by 

w , x s _ 7r f (xi, . . . ,a: f ) 7r t (x t ) Lt-i {x t ,x t -i) 

The resulting particles are now weighted samples from 7r t . Consequently from ([]]), 
under the SMC sampler framework, one may work directly with the marginal dis- 
tributions 7i t (x t ) such that w t (xi, . . . , x t ) = Wt(x t -i,x t ). While the choice of the 
backward kernels L t _\ is essentially arbitrary, their specification can strongly affect 



the performance of the algorithm. See Del Moral et al. (2006) for detailed discussion 



2.2 Incorporating partial rejection control 

It is well known that the performance of SMC methods are strongly dependent on 
the mutation kernel (ICornebise et al. 2008|) . If M t is poorly chosen, such that it 



does not place particles in regions of the support of ix t with high density, then many 
importance sampling weights will be close to zero. This leads to sample degeneracy, 
as a few well located particles with large weights dominate the particle population, 
resulting in large variance for estimates made using these samples. 



Liu (2001) (see also lLiu et al. 19981) introduced the partial rejection control strat- 
egy to overcome particle degeneracy in a sequential importance sampling setting. 
Under this mechanism, when the weight of a particle at distribution ir t falls below 
a finite threshold, c t > 0, the particle is probabilistically discarded. It is replaced 
with a particle drawn from the previous distribution ir t _i which is then mutated to 
%f This new particle's weight is then compared to the threshold, with this process 
repeating until a particle is accepted. This approach is termed partial rejection, as 
the replacement particle is drawn from ir t _i, not ~K\ (ILiu et al. 19980 . 



Under the SMC sampler framework we modify this approach and incorporate 
the partial rejection mechanism directly within the mutation kernel. Hence at time 
t — 1, the particle xt-i is moved via the mutation kernel M t (xt-i,Xt) and weighted 
according to (fTj) . This particle is accepted with probability p, determined by the 
particle's weight and the weight threshold c t . If rejected, a new particle is obtained 
via the mutation kernel M t , until a particle is accepted. 

For the sequence of distributions ir t , t = 1, . . . , T, the mutation and backward ker- 
nels Mt and L t -i, a sequence of weight thresholds q, and PRC normalizing constants 
r(c t ,x t -i) (defined below), the SMC sampler PRC algorithm is given by: 
SMC sampler PRC algorithm 

Initialization: Set t = 1. 

For i — 1, . . . , N, sample x± ~ 7Ti (ar) , and set weights Wi(xi ) = 4. 



Resample: Normalize the weights Y^i^t{ x t ) = 1- If Ei Wt resam- 



(b) With probability 1 — = 1 — min{l, Wt{x\ )/c 4 }, reject x[ and go to (a). 

(c) Otherwise, accept and set W t (x[^) = W t (x^)r(c t , x^_^)/p^. 

(d) Increment % = % + 1. If % < N, go to (a). 

(e) If t < T go to Resample. 

The above algorithm without the mutation and correction steps (b) and (c) is 
equivalent to the standard SMC sampler algorithm (IDel Moral et al. 20060 . In the 
resample stage, the degeneracy of the particle approximation is quantified through the 



pie N particles with respect to W^xf 1 ) and set Wjt(x^) 




1,...,JV. 



Mutation and correction: Set t = t + 1 and % — 1: 

(a) Sample x[ 1 ' ~ M^x^l^Xt) and set weight for x^ to 
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usual estimate of the effective sample size, 1 < \^2 i W t (x[ l ~ > ) 2 ]~ 1 < N (ILiu and Chen 19981) 
We discuss the choice of the thresholds, c t , in later Sections. 

The addition of a rejection step at each time t effectively modifies the mutation 
kernel M t . We denote the resulting kernel by M£, where 



M* (x t -i,x t ) = r(ct, x t -i) 1 min 



Ct 



• (2) 



The quantity r(c t ,x t -i) denotes the normalizing constant for particle Xt-i, given by 

r(c t , x t -i) = J min jl, W t -\{x t -i) — | M t (x t -i, x t ) dx t . (3) 

Note that < r(ct,Xt-i) < 1 if (w.l.o.g.) the mutation kernel M t is normalized, so 
that J M t (xt-i,Xt)dxt = 1, and if the PRC threshold < c t < oo is finite. Thus the 
SMC sampler PRC algorithm can be considered as an SMC sampler algorithm with 
the mutation kernel M* (xt-i,Xt), and the correction weight 

rrr, \ T/f/ ( \ K t (x t ) L t _l (x t , X t -l) . , 

2.3 Estimation of the normalizing constant 

As the normalizing constant r(ct, Xt-i) in the weight calculation (jl]) in general depends 
on x t -\i it must be evaluated. Where no analytic solution can be found, approximat- 
ing (jBD may be achieved by, for example, quadrature methods if the sample space E 
is relatively low dimensional or Monte Carlo methods if E is high dimensional. For 
example, for j = 1, . . . , m independent samples x* t ^ sampled from M t {x t -i^x t ) 



1 ™ w \ X t -i,X t ' 

r(ct,x t -i) w — > min { 1, W t _ 1 (a;t_ 1 )- 
i=i { 

An alternative, computationally more efficient approach is to select kernels M t and 
L t _i such that r(c u x t -i) = r(c t ) will be constant for all particles Xt-i- In this case, 



the value of r(c t ) may be absorbed into the proportionality constant of the weights, 
and safely ignored. Equation suggests that this can be achieved if M t (x t -i,x t ), 
W t -i(x t -i) and w(xt-i,Xt) are independent of Xt-\- 

Specifying mutation kernels M t such that M t (x t -i, x t ) = M t {x t ) amounts to choos- 
ing a global kernel which is the same for all particles Xt-i- This is common in practice 
(e.g. IWest 19 93). The particle dependent weight W t -i(x t -i) can be set to 1/N for all 



particles following a resampling (or preselection) step. Finally, consider for a moment 
the backward kernel of the form 

rapt ( v _ 7r t -i(xt-i)M t (xt-i,Xt) . s 

J 7rt-i{xt-i)Mt{xt-i,Xt)dx t -i 

This backward kernel is an approximation of the optimal backward kernel, in the sense 
of the choice of L t -\ that minimizes the variance of the importance sampling weights 
(IDel Moral et al. 2006|) . Under the backward kernel (jSj), the incremental weight can 



be approximated by 

w t (xt-uxt) = irt(x t )/ I n t -i(xt-i)Mt(xt-i,xt)dx t - 



N 



(xtyYsW^ixfljMtixfl^Xt). 



1=1 

Under a global mutation kernel M t (xt), and following a resampling step, then the in- 
cremental weight under this backward kernel reduces to Wt(x t -i,x t ) = ir t (x t ) / Mt(x t ) , 
which is independent of x t -\- Thus, the weight calculation in (j3J) becomes 



W t (x t ) oc n t (x t )/ 



mm 



w(x t -i,xt) 



■K t {x t ) / M t (x t ) ifminjl, Nct 
Net otherwise. 

Note that under this setting, the SMC sampler PRC algorithm can be considered as 
a sequence of importance sampling strategies with partial rejection control. 



3 SMC Sampler PRC algorithm analysis 

In this section we study theoretical properties of the SMC sampler PRC algorithm. 
We firstly bound the variance of the importance weights, and then present a central 
limit theorem for the sampler with a recursive expression for the asymptotic variance. 
Finally, via a connection with an existing algorithm, we establish a condition for which 
the number of rejection steps under the PRC mechanism is almost surely finite. 

3.1 Variance of the incremental weights 

We begin this section by establishing a bound on the variance of the importance 
weights of the SMC sampler PRC algorithm. 

Theorem 3.1 Let W t {x t ) denote the importance sampling weight at time t from a 
standard SMC sampler with mutation kernel M t , and let (x t ) denote the equiv- 
alent weight following a partial rejection control step under the SMC sampler PRC 
algorithm, with resulting mutation kernel M£ . Then 

Yar M * [W; (x t )} < Var Mt [W t (x t )} . 

Proof: See Appendix A.l. 

Hence, applying partial rejection control within the SMC sampler framework will 
not worsen, and may improve the variance of the importance weights, by reducing 
the x 2 -distance between the sampling and target distributions at each stage, t. 

In the case where min jl, Wt ^ j = 1 for all x t , which is achieved when c t < 
mf Xt {W t (x t )}, then from (jHJ) we have r(c t ,x t -i) = 1 for all x t -i- From (j2J), this 
results in M^(x t -i,x t ) = M t (x t -i,x t ) and hence Yar M * [W^ (x t )} = Yar Mt [W t (x t )]. 
That is, the SMC sampler PRC algorithm reduces to the standard SMC sampler 
when ct < inf Xt {Wt (xt)}, and in this case, the variance of the importance weights is 
maximised. When W t (x t ) G [0, oo) this is realized for c t = where we define 0/0 := 1. 
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3.2 A central limit theorem 

Central Limit Theorems (CLTs) for SMC and particle filtering algorithms have been 

derived in various literatures fIDel Moral 20041 IKunsch 20051 |Chopin 2004] IDel Moral et al. 2006 



IJohansen and Doucet 20081) . They are based on the observation that an SMC algo- 
rithm introduces local errors (fluctuations) as a result of the approximations intro- 
duced by sampling numerically from the transitions. Hence, at each stage t, one can 
decompose the error between the target distribution ir t and the iV-particle approx- 
imation . This turns out to be a sum of the local sampling fluctuations at each 
discrete time in the past, propagated forwards in time to t. 

In the setting of the SMC sampler algorithm, the existence of a CLT is established 



by Del Moral (2004) Explicitly, under the assumption of multinomial resampling at 



each stage of the algorithm, and the integrability conditions given in Chopin (2004) 



[Theorem 1] and Del Moral (2004) [Section 9.4, pp. 300-306], then for a suitable con- 



tinuous and bounded test function if G Cb (E) we have 

N 1 ' 2 (E^) - E^y?)) - JV(0, VsmcM) ( 6 ) 



as N — > oo, for each t = 1, . . . , T. Del Moral et al. (2006) obtain a recursive expres 



sion for the asymptotic variance VsMC,t(f) as an explicit function of the backward 
kernels L t -\ and the sequence of distributions on path space, 7r t - 



Following Del Moral et al. (2006) , we obtain an analogous result for the SMC 
sampler PRC algorithm. Under the same assumptions as the above, we have the 
CLT ([6]) with asymptotic variance given by 

VsMC-PRC,t (tp) = J h^^j (J V? Oi) tti (ztki) dx t - E nt (p)) 2 dxx 



t-i 



k=2 
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with h. 



-i 



and 1 1 = 1. 



(c fc ,x fc _i) mm jl, iv^.^.OM.^^c, J 
The contribution of the PRC stage to the asymptotic variance is encapsulated in 
the Ik terms. Under the standard SMC sampler algorithm we have Ck < inf Xt {Wf (xt)} 
so that Ik — 1 for all k — 1, . . . , t. In this setting, ((Tj) reduces to the asymptotic 



variance expression obtained by Del Moral et al. (2006) 



3.3 Connections to an existing SMC algorithm 

In rare event applications there is a high probability of generating particles with ex- 



actly zero weights. The AliveSMC algorithm (Le Gland and Oudjane 2004 Le Gland and Oudjane 2005) 
was developed to ensure that a particle population of a desired size persists at each it- 
eration of a standard SMC algorithm (see lDel Moral et al. 200i~| IJohansen et al. 20061 



for related methods). In this setting, the number of particles at each time t is con- 
sidered as a random variable N Ct . That is, N Ct is the number of particles required 
to generate exactly N non-zero weighted particles. In this Section we reinterpret the 
SMC sampler PRC algorithm in terms of the AliveSMC algorithm. As a consequence, 
in Section 13.41 we are able to establish a condition under which the PRC resampling 
stage will require a finite number of rejection attempts. 



In Le Gland and Oudjane (2005), at iteration t, a fitness function is applied to 
select particles satisfying a desired criteria. Those particles not satisfying the criteria 
receive a zero weight. In an SMC sampler PRC setting, the fitness function can be 
interpreted as selecting those particles with a weight that is immediately accepted 
under the PRC acceptance probability. As such, we may rewrite the SMC sampler 
PRC algorithm with a modified mutation and correction step: 
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Reinterpreted SMC sampler PRC/ AliveSMC algorithm 



Mutation and correction: Set t — t + 1 and j — 1. For i = 1, . . . , N ( 



C1 ■ 



W t {xf) oc 



(a) Sample xf ~ M t (x® u x t ) and calculate W = W^x^) Mx * u^fl^'w • 

(b) Set weight for x[ as 

Wr(ct, x t _i)/p^ with probability — min{l, W/ct} 
otherwise. 

(c) If Wt(x^) 7^ 0, increment j = j + 1. 
If t < T go to Resample. 

Note that j — 1, . . . , N indexes the particles x^_i at time t — 1 such that particle 
mutations from x^}\ generate a non-zero weight exactly once. Also, under the fitness 
function, particle xf 1 has a probability 1 — = 1 — min{l, W/ct} of being exactly 
zero. Given that it is possible to express the SMC sampler PRC algorithm within the 



AliveSMC framework, we may adapt the results of Le Gland and Oudjane (2004) and 



Le Gland and Oudjane (2005), to obtain a condition under which the SMC sampler 
PRC algorithm is guaranteed to require a finite number of attempts, N Ct < oo, to 
obtain exactly TV non-zero weighted particles. 

3.4 Analysis of the number of rejection attempts 



Following Le Gland and Oudjane (2004) and Le Gland and Oudjane (2005), we de 



fine the random variable N Ct as 

N* 



N ct = inf { N* ^ 1 : V W^(x t ) > N sup W t (x t ) 

xt&E 



t=l 



Le Gland and Oudjane (2005) proved for the AliveSMC algorithm that the random 
number of particles N Ct is almost surely finite with N Ct ^ N, under the condition 

St) / Tr t -i(x t -i)Mt(xt-i,xt)d 
f ir t -i(xt-i)Mt(xt-i,xt)dxt- 
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that (TTt-iM t ,W t ) = 1 Wt{ i 1 ^ {x \^f t{xt -\ x ; )dxt ~ ldxt > 0. A sufficient condition for 

» ■* ± *' k ' 7T+—1 If- 1 Alt (Xt — 1 .xt )dxt--\ 



this to hold is E Mt [W t (x t ) \ x t -\ — x] > 0, for all x G E. Thus for the SMC sampler 
PRC algorithm, we have the following theorem: 

Theorem 3.2 Under the SMC sampler PRC algorithm, the number of rejection at- 
tempts at each stage of the algorithm, N Ct N , is almost surely finite if c t < oo. 

Proof: See Appendix A. 2. 

Corollary 5.1 The following convergence in probability holds, with a rate of 1/y/N: 



N (7r t _iM tj W t ) 



< oo. 



See Le Gland and Oudjane (2005) for further details. Hence, the SMC sampler PRC 
algorithm possesses an almost surely finite number of rejection attempts if the PRC 
threshold c t is finite, with the above rate of convergence. 



4 Approximate Bayesian computation 

With the aim of posterior simulation from tx{x\T>) oc ti{T>\x)ti{x) for parameters x 
and observed data T>, "likelihood-free," approximate Bayesian computation (ABC) 
methods are often utilised when the likelihood function, ir(T>\x), is computationally 
intractable or when its evaluation is computationally prohibitive. ABC methods can 
be based on rejection sampling (ITavare et al. 1997j IBeaumont et al. 20021) . Markov 
chain Monte Carlo (|Marjoram et al. 2003} IBortot et al. 20071 IRatmann et al. 20091) 



and SMC-type samplers fISisson et al. 20071 IToni et al. 20091 IBeaumont et al. 20091 
IDel Moral et al. 20081) . While currently among the most efficient ABC methods, the 
underlying practical issue with SMC-type algorithms is in avoiding sample degeneracy 
through extreme numbers of particles with low or exactly zero weights. In this Section, 
we will demonstrate that the SMC sampler PRC algorithm applied within the ABC 
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framework can achieve significant performance gains and greater modelling flexibility 
over existing SMC-type ABC methods. 

The underlying approach of ABC methods is to augment the (intractable) pos- 
terior to n(x, T>'\T>) oc tt(T>\T>', x)tt('D , \x)'k(x), where the auxiliary parameter is an 
artificial data set distributed according to the model V ~ tt(-\x). An approximation 
of the target posterior tt(x\V) is then given by 



The weighting distribution n(V\V, x) takes high density in regions where the datasets 
T> and T>' are similar, and low density otherwise. Comparison of the datasets is usually 
achieved through low-dimensional summary statistics T(-), so that, for example 



for some small tolerance value e > and distance measure p. If T(-) are sufficient 
statistics, and e — > so that tt(V\V,x) reduces to a point mass at T(T>) = T{T>') 
then n abc { x I — t^(x\T>) is recovered exactly, otherwise the ABC approximation to 
7r(x\D) is of the form ([8]), with greater accuracy for smaller e. The computational 
overhead of all ABC samplers increases as e decreases, producing a trade off be- 
tween computation and accuracy. ABC methods either sample from the joint density 
tt(x, T>'\T>) by arranging to cancel out the intractable likelihood in a weight or accep- 
tance probability, or sample from ttabc{x\D) directly via Monte Carlo integration 



where T>[, . . . , T>' s ~ %(T>'\x) are draws from the likelihood given x. Almost all current 
ABC methods have the weighting density ([9]) written directly into the algorithm. 

We apply the SMC sampler PRC algorithm in the ABC framework as follows: The 
target 7i t (x t ) = ^ABCt^t^) is given by ([8]), with the weighting function n t (V\V',x) 





(9) 




(10) 



s=l 
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parameterized by a different scaling parameter e t for each t, where oo = e\ > . . . > ex, 
produces increasing accuracy at each step, t. The et sequence and its length, T, may 
be determined a priori or dynamically. Evaluation of 7r t (xt) is defined by ([TO]) through 
S Monte Carlo draws. Given the high computational overheads of ABC methods, we 
avoid evaluating the PRC normalizing constant (as demonstrated in Section 12.3ft . 
through a global mutation kernel M t (x t ), the backward kernel L°^\ (c.f. and 
enforced resampling. 
SMC sampler PRC-ABC algorithm 

Initialization: Set t — 1. 

For i — 1, . . . , N, sample Xi ~ fi{x), and set weights Wt(x^) = n abc,i{ x i^ '|^) / l^( x i^) ■ 

Resample: Resample N particles with respect to W t (x{ ) and set W t (Xf ) = j^, 
i = l,...,N. 

Mutation and correction: Set t — t + 1 and % — 1: 

(a) Sample x^ ~ M t [xt) and set weight for xf 1 to 

W t {xf) = 7i ABC ,{xf\V)/M t {xf). 

(b) With probability 1 — = 1 — min{l, Wt(xf j/c*}, reject aj^ and go to (a). 

(c) Otherwise, accept and set Wt(x\ ) = Wt(x®)/P® • 

(d) Increment % — % + 1. If z < iV, go to (a). 

(e) If t < T then go to Resample. 

The density n(x) is an initial sampling distribution, from which direct sampling 
is available. As with the tolerance e t , the PRC thresholds c t may also be deter- 
mined dynamically (see below for an illustration). Note that as the resampled par- 
ticles in the Resample step play no subsequent part in the sampler, in practice this 
step can be omitted. The path of each particle (x± \ . . . , x^ ) £ E T can be recon- 
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structed post-simulation, if required, by resampling the recorded marginal populations 

{W t (xf),xf). 

The above algorithm has a number of benefits over existing SMC-type ABC sam- 
plers fISisson et al. 20071 IToni et al. 20091 IBeaumont et al. 20091 IDel Moral et al. 2008j) 



Firstly, the weighting density tt(T>\V, x) can take any form - we suggest any smooth- 



ing kernel, following Blum (2009) Existing samplers in the literature are restricted 



to the uniform function (Q. Secondly, there is complete control over the PRC thresh- 



old, q, unlike Sisson et al. (2007) who impose a specific value. Thirdly, in estimating 
' K ABc( x t\'D)i a s long as the backward kernel is used, any number S > 1 of Monte 
Carlo draws can be used in (jTOj) . Existing samplers only use S = 1, and so there is less 
control over the variability of the weights. Finally, providing that the computation 
required to estimate the PRC normalizing constants, r(c t , x t -i), is acceptable, a form 
of the SMC sampler PRC- ABC sampler may be constructed which uses arbitrary mu- 
tation and backward kernels, allowing the user to select the most appropriate tools 
for a given problem. 



4.1 Simulation study 

We now demonstrate the superior performance of the SMC sampler PRC-ABC algo- 
rithm through a controlled study. Specifically, we specify the true posterior n(x\V) 
as iV(0, 1) by defining the likelihood and prior as T> ~ N(x, 1) and ir(x) oc 1, with 
a single observed datum, T> — 0. For this model, a sufficient statistic is T{T>) = T>. 
From (El), for the uniform weighting density with p(a, b) = \a — b\ and e = e u , or 
for tt(T>' I T>, x) = iV(2),ep, then tiabc{ x \T^) ma Y be obtained in closed form as 

7r ABC (x\V) oc $(£ " ~ X) ~ H ~ €u - X) or tt abc (x\V) = N(0, 1 + e 2 q ) 
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respectively, where $(•) denotes the standard Gaussian CDF. In both cases ttabc( x \'D) —> 
N(0, 1) as e — > 0. In order to directly compare the two approximate posteriors we 
impose equal variances on the two weighting functions, so that e g = V3e u 

We adopt the following sampler specifications: A particle population of size N = 
1000 was drawn from the initial sampling distribution fi(x) ~ U(— 5,5), and the se- 
quence of distributions, ir±, . . . , 7r 10 , is defined by {e t } = {oo, 10, 5, 2, 1, 0.5, 0.2, 0.1, 0.05, 0.05}, 
on the e g scale. The mutation kernel M t (xt) = Ylf=i ^t-i l l ) { x t\ x 't-ii r2 ) * s taken as 
a Normal kernel density estimate of 7T t _i(x t -i), with r 2 = 1 and where if)(x\fi,o~ 2 ) 
denotes the PDF of a N(fi, a 2 ) distribution evaluated at x. We initially use S = 1 
Monte Carlo draws to approximate •kabc( x \'D) (c.f. [TO]) . 

Figure [TJ examines the effect of PRC on the effective sample size (ESS), the 
variance of the importance weights and the mean number of rejections per par- 
ticle. The PRC threshold was determined dynamically at each iteration as c t = 
Q(W^~(x t ),q), the g-th quantile of the non-zero weights at time t (obtained by mu- 
tating all Xt-i particles under M t before implementing the PRC stage), for q = 
0,0.5,0.75,0.85,0.9,0.95,0.99,0.995,0.999. Results are shown using the Gaussian 
(left plots) and uniform (right plots) weighting density, based on 250 sampler repli- 
cations. Note that the PRC threshold with q = approximately corresponds to a 
standard SMC sampler ( "No-PRC" ) only for the Gaussian weighting function, as the 
uniform weighting function permits exactly zero importance weights. Setting q = 
for the uniform weighting density corresponds to existing SMC-type ABC samplers 
( ISisson et al. 20071 iTbm et al. 20091 Beaumont et al. 2009|) . 

For both weighting densities, the effective sample size increases as c t increases, 
and the variance of the importance weights decreases. Naturally, the higher the PRC 
threshold, the more rejections occur, quantifying the extra computation required for 
the gains in sampler performance. However, there is a notable difference in the 
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(a) (b) 




No-PRC 0.5 0.75 0.85 0.9 0.95 0.99 0.995 0.999 No-PRC 0.5 0.75 0.85 0.9 0.95 0.99 0.995 0.999 

PRC threshold quantile PRC threshold quantile 



Figure 1: Effective sample size (a,b), variance of normalized importance weights (c,d) 
and mean number of rejections per particle (e,f) as functions of PRC threshold ct- PRC 
threshold is defined dynamically as a quantile of the non-zero weights at time t (x-axis). Left 
plots (a,c,e) and right plots (b,d,f) are obtained under the Gaussian and uniform weighting 
densities i:(V\V ,x) respectively. Boxplots are based on 250 sampler replications. 

transition from poor (no PRC) to improved (under PRC) performance between the 
two different weighting densities. This occurs as the uniform weighting density only 
permits 0/1 weights, compared to the smoother scale under the Gaussian. As a result, 
the uniform weighting density (which is the only choice under existing ABC samplers) 
has a fixed, albeit strong, performance gain over not implementing PRC, but at a 
very high computational cost (Figure (Uf). Comparison with panels (a,c,e) suggests 
that considerable computational gains can be achieved with alternative weighting 
functions, without sacrificing sampler performance. This is easily permitted under 
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the SMC sampler PRC framework. 



(a) 



(b) 



10 50 100 500 1000 2000 

Number ol Monle Carlo draws 

(C) 



10 50 100 500 1000 2000 

Number ol Monle Carlo draws 



50 100 500 1000 2000 

Number of Monle Carlo draws 

(d) 



10 50 100 500 1000 2000 

Number of Monle Carlo draws 



Figure 2: The effect of the number of Monte Carlo draws, S, on sampler performance. 
Panels show (a) effective sample size (ESS), (b) variance of the importance weights, (c) the 
mean number of rejection attempts, and (d) estimates of the posterior variance (true value 
~ 1), as a function of the number of Monte Carlo draws in the estimation of kabc{ x \D)- 



When using the L°^_ x backward kernel (jSJ), any number S > 1 of Monte Carlo 
draws may be used to approximate ttabc{ x \'D) via ffTUl) . While S — 1 is near universal 
under existing ABC algorithms, one would expect to realize less variable importance 
weights for S > 1. Figure (j2j) illustrates the effect of increasing S, under the Gaussian 
density function, based on the PRC threshold c t = Q{W^{xt), 0.95). An increase in 
the effective sample size (panel a) is reflected by the reduction in the variability of the 
importance weights (panel b), as is the variability in the estimates of the posterior 
variance (panel d). This in turn results in lower numbers of rejections at the PRC 
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stage (panel c). Of course, these performance gains are again balanced by the strong 
increases in computation required for S > 1. It would appear that unless the data- 
generation procedure T>' s ~ 7r(T>\x) is computationally inexpensive, 1 < S < 10 
would seem to be the most useful choice in practice. Regardless, the greatest gains in 
sampler performance under the SMC sampler PRC algorithm are achieved for S = 1. 

5 A stochastic claims reserving analysis 

We present an analysis of an important and popular class of statistical models in actu- 
arial science using stochastic claims reserving. We consider a time series formulation 
of the distribution-free chain ladder model (IMack 1993( IGisler and Wiithrich 2008| 
IPeters et al. 20081) . For a claim on an insurance company for an accident in year i, 
Cij denotes the cumulative claim in subsequent years j > i. Cumulative claims can 
refer to payments, claims incurred and other expenses. At time J, we have obser- 
vations T>i = {Cij; i + j < I}, and for reserving against future claims we wish to 
predict = {Cy; i + j > I, i < I}. One such dataset is illustrated in Tabled] 

Under a time series formulation, cumulative claims in different accident years 
i are independent and satisfy, for j = 0, 1, 

Cjj'+l = fjCi,j + °j \J Ci,j e i,j+1> (11) 

where f = (/ , • • • , fi-i) and cr = (cr , . . . , 07_i) are respectively the chain ladder fac- 
tors and standard deviations, and the residuals are i.i.d. with mean and variance 
1. The model is constrained such that P(Cij > 0\{Ck,o} 3 k = 1 ,f,<r) = 1 for all i,j (see 
IPeters et al. 20081) . If distributional assumptions are made on the residuals e^j (e.g. 
lYao 20081) . the posterior distribution can be made computationally tractable. How- 
ever, a primary intention of this model is to work with distribution-free assumptions 
on the residuals, and therefore on the cumulative claims. Within this distribution-free 
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context one wishes to quantify popular risk metrics such as value-at-risk and expected- 
shortfall to be calculated for the predicted claims distribution, both of which are 
highly relevant to regulatory reporting. Alternative approaches, based on credibility 
results, can relax such distributional assumptions, but can only provide statements 
on posterior first and second moments in limiting cases (IGisler and Wuthrich 20081) . 

Previously, actuaries have proceeded by predicting claims via a deterministic 
model known as the classical chain ladder algorithm. This approach predicts un- 
observed future cumulative claims by the recursion Cij-i = Cij-i, and for j > I — i 

C id = where gg? = jgff , (12) 

and where, in the time series formulation, the variances are estimated by 



/ \ 2 



I , 1 I C id Ij 



j 

-2(CL) 



See Mack (1993) for an estimator of o i I _ l . As this algorithm is deterministic there is 
strong interest in stochastic chain ladder models, which naturally allow the quantifica- 
tion of uncertainty, such as the mean square error of prediction. In the claims reserving 
setting the most popular stochastic models are those with estimators which recover 
the classical chain ladder estimators. We consider one such Bayesian stochastic model 
which has the property that as the diffusivity of the priors 7r(f , cr) tends to infinity 
f(MMSE) f (CL) w h ere MMSE denotes the posterior mean (IGisler and Wuthrich 200£ 
Hence by (1121) the posterior mean E[Ci : j\Vj] = C^j recovers the classical estimators, 
thereby justifying the classical model. We sample from the intractable posterior 
7TABc(f , cl'C/) using the SMC sampler PRC- ABC algorithm. 

5.1 Analysis and results 

This model is interesting as the intractability of the likelihood directly impacts the 
ability to generate synthetic data sets, T>\ 1 from the model. That is, if the distri- 
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butional form of the residuals were known, data-generation from the model would 
be trivial. To retain a distribution-free setting we alternatively utilise a conditional 
bootstrap approach (IPeters et al. 20081) . Conditional upon proposed parameters f 
and er, the residuals eij|f, <r are iteratively obtained by inversion of (fTTT) . Then, 
by independently drawing resampled residuals from the empirical conditional resid- 
ual distribution, a bootstrap sample of the cumulative claims D\ is then available 
through recursion on (fTTj) . 

In analyzing the real claims reserving data in Table Q] we specify independent 
priors fj ~ Gamma(aj, (3j) with mean aj(3j = fj CL ^ and aj ~ IGamma(a.,-, bj) with 
mean bj/ (aj — 1) = aj CL ^ for j — 0, . . . , I — 1, each with large variance. For summary 
statistics we adopt T{V') = (V'j, //(e), s'(e)) where //(e) and s'(e) denote the sample 
mean and standard deviation of the conditionally resampled residuals e^|f,cr. The 
observed summary statistics are given by T(T>) = (Vj, 0,1) following the zero mean 
and unit variance assumptions on the true residuals. 

We implement the SMC sampler PRC-ABC algorithm with uniform weighting 
density © and p(T(Vr), T{V'j)) = [(T(V I )—T{V' I )) J T I ~ 1 {T(V I ) —TiV'j))] 1 / 2 defined 



as Mahalanobis distance, where the covariance £ is estimated following Peters et al. (2008) 



We use iV = 5000 particles, PRC threshold c t = Q(W^(x t ), 0) and a deterministic 
distribution schedule {e t } = {oo, 10, ... , 0.00001} with T = 22. The mutation kernel 
M t (x t ) = YliLi W / t-''i( x i-i)G'amma(a(x|^ 1 ), b(xfl 1 )) is a mixture of gamma densities, 
with mean a(x^ 1 )6(x{^ 1 ) = xf]_i and large variance. 

Table [2] presents a comparison of the parameter estimates f and <x, and predicted 
cumulative claims, Cjj under classical and Bayesian models. Given the uninformative 
priors, the posterior mean estimates and resulting predicted claims agree well with 
those obtained under the classical model. This provides some validation for the 
deterministic classical model estimates under the Bayesian stochastic interpretation. 
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Perhaps more usefully for inference, the full posterior / Kabc{^i (J '\T^i) is available. 
Figure [3] illustrates how the estimated marginal densities of the first chain ladder 
factor, 7iABC,t(fo\'E > i), and the associated standard deviation, iiABC,t(o'o\'Di), evolve as 
e t decreases. The precision of the densities clearly improves, as decreasing e t imposes 
stricter restrictions on the permissible deviations of the ABC approximate posterior 
7TABc(f > 47 1 ^) from the target posterior 7r(f , cr\T>). A full predictive analysis may now 
follow, including upper and lower credible bounds on predicted future claims. 



6 Discussion 

When used in challenging settings, sequential Monte Carlo samplers often suffer from 
severe particle degeneracy. In this article we have provided a practical approach 
to tackling this problem, by incorporating the partial rejection control mechanism 



of Liu (2001) directly into the mutation kernel of the SMC sampler. The resulting 
sampler will not worsen, and can improve the variance of the importance weights, 
sometimes substantially so. By establishing clear relationships with existing sam- 



plers (IDel Moral et al. 20061 |Le Gland and Oudj ane 2004|), many theoretical proper 



ties may be extended to the SMC sampler PRC algorithm, including a central limit 
theorem, and a proof of an almost sure finite number of PRC rejection attempts. 

There is much opportunity for the specification of the sequence of PRC thresholds 
to be further automated, if desired. For example, by dynamically determining q > 
if the effective sample size of the particle population at time t — 1 falls too low, 
and conversely allowing c t < c t -\ if the level of PRC resampling is too high, in order 
to reduce computational overheads. 

As the SMC sampler PRC algorithm allows practical inference in challenging sit- 
uations in which particle weights are highly variable, we anticipate that a primary 
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application of the sampler will be within the rapidly developing "likelihood-free" ap- 
proximate Bayesian computation framework. The presented sampler is more flexible 
and efficient than existing SMC-type ABC samplers, allowing a previously unavailable 
degree of control over the computation utilised for a given analysis. Perhaps more 
importantly, the extra flexibility achieved by allowing arbitrary weighting densities 
(unlike existing ABC samplers) enables the analysis of improved models within the 
ABC framework, in line with recent non-parametric interpretations (IBlum 2009|) . 
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Appendix 

A.l: Proof of Theorem SH] 



The proof follows the arguments presented by Liu et al. (1998) In particular we 
study the SMC sampler PRC algorithm in terms of x 2 distance between the sampling 
distribution and the target distribution at stage t. Let 

W t (x t ) oc Wt-\(x t -i) -, r . _ , r and W t (x t ) oc Wt-i[x t -\) 



7Tt_i {Xt-l) M t (x t -l,X t ) TTt^l (Xf_l)M t * (x t -l,X t ) 

Recall that the normalising constant for the mutation kernel Ml at time t is 



J min |l, | M t (x t -i, x t ) dx t = ^E Mt [min {c t , W t (x t )}] ■ 



The variance of the importance weight at time t from a standard SMC sampler, with 
respect to M t , is given by 

Yar Mt [W t (x t )] = J [W t (x t )] 2 M t (x t ^ x t ) dx t - 

and similarly, the variance of the equivalent importance weight at time t following a 
PRC step under the SMC sampler PRC algorithm, with respect to M t *, is given by 

Var Mt * [W*\x t )\ = J [W;{x t )fNi; (x t -x, x t ) dx t - /i 2 , 

where 

/2 = E M ^[W t *(x t )} =E Mt [W t (x t )}. 
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We also have that 

2 



f[W t *( Xt )} 2 M*(x t ^x t )dx t = J 



Tr{x t )L t -i(xt,xt-i) 



M* (x t -l,Xt)dXt 



7T t _l(x t _l)M*(x t _l,Xi) 

r ^ X *-V J —TILSIT ^TA^JWi^^ t[Xt ~ u Xt)dXt 



W^xt.r) ^(x t )L 2 t _ 1 (x t . 1 ) 

r(ct,x t -i) f raa,y.{Wl(x t ),c t W t (x t )}M t (x t _i,x t )dx t 



= r(ct,Xt-x)E Mt {m&x{W t (x t ),Ct}W t (x t )} 

= j-E Mt [min fa, W t {x t )}\ E Mt [msx{W t (x t ), c t }W t {x t )} 

< }E Mt [min {a, W t (x t )} mnx{W t (x t ), c t }W t {x t )} 

= ±E Mt [c t W 2 (x t )} = E Mt [W 2 (x t )} . 

The above inequality holds since the random variables min {q, Wt(xt)} and m&x{Wt(xt), Ct}Wt 
are positively correlated (see ILiu 2001]) , and so 



E Mt [min {c u W t (x t )} max{W t (x t ) , c t }W t (x t )} 

-E Mt [min{c t ,W t (^)}] E M( [max{W t (x t ),Ct}W t (x t )} > 0. 

Hence 

Var M ;\W t *(x t )] < E Mt [W 2 (x t )] - fi 2 = ¥ar Mt [W t (x t )}. 

□ 



A. 2: Proof of Theorem SH 

To satisfy the condition E Mt [W t (x t ) \ x t -i = x] > 0, Vx G E, for the SMC sampler 
PRC algorithm, we require 

/ Wt-i(x)w t (x,x t )r(c u x) 

J E 

For a particle xt-i = x, the proposed state Xt can take values in a support which 
can be split into two regions, A(x) and A c (x), such that A(x) U A c (x) = E. These 
respectively correspond to when min jl, j = \ ^ e particle acceptance 

probability under PRC is 1) and when min jl, ^'-'M^t 1 ' 1 ') j — ty t _ 1 (x)w t (x, x t ) /c t 
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min < 1 



M t (x,x t )dx t >0. (13) 



(i.e. the particle may be rejected under PRC). Note that in the extreme cases, 
Q < Wt~i(x)w t (x, x t ) =>- A c (x) = reduces to the SMC sampler algorithm, and 
c t > sup Xt Wt-i(x)wt(x, x t ) =>• A(x) = 0. More generally, f|T3|) may be expanded as 

v(ct x)Wt-i(x) f f 

— / ix t (x t ) L t -i (x t , x) dx t + c t r(c t , x) I M t (x, x t )dx t > 

7Tt-l {%) JA(x) JA c {x) 

which is always greater than zero for finite c f < oo as < r(c t ,x) < 1. 

□ 



Accident 
Year, i 





1 


2 


Development Year, j 
3 4 5 


6 


7 


8 


9 





594.6975 


372.1236 


89.5717 


20.7760 


20.6704 


6.2124 


6.5813 


1.4850 


1.1130 


1.5813 


1 


634.6756 


324.6406 


72.3222 


15.1797 


6.7824 


3.6603 


5.2752 


1.1186 


1.1646 




2 


626.9090 


297.6223 


84.7053 


26.2768 


15.2703 


6.5444 


5.3545 


0.8924 






3 


586.3015 


268.3224 


72.2532 


19.0653 


13.2976 


8.8340 


4.3329 








4 


577.8885 


274.5229 


65.3894 


27.3395 


23.0288 


10.5224 










5 


618.4793 


282.8338 


57.2765 


24.4899 


10.4957 












6 


560.0184 


289.3207 


56.3114 


22.5517 














7 


528.8066 


244.0103 


52.8043 






to be predicted Yij 






8 


529.0793 


235.7936 


















9 


567.5568 





















Table 1: A claims development triangle. Upper triangle denotes observed annual claims 
Yij from which Cjj = YjL=o ^ ms ¥ be obtained; lower triangle denotes annual Y{ j 
and cumulative claims Cij £ T>j to be predicted. Data are real insurance figures in units of 
$10,000 (cf. Wiithrich and Merz, 2008). The triangle assumes that the number of accident 
years is equal to the number of observed development periods. 
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Figure 3: Evolution of the marginal posterior density estimates of the chain ladder factor 
ir t (fo\Vj) (left) and the associated standard deviation 7Tt(ao\T>i) (right) as a function of e<. 
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9, 862, 404 


9, 877, 198 


90, 283 


(f< CL )) 


5 












10, 005, 044 


10, 056, 528 


10, 067, 393 


10,077,931 


10, 092, 247 


156, 494 


1 Y\ 3£i J i 














10,019,212 


10,074,318 


10,087,415 


10, 096, 493 


10, 111,638 


175,886 




6 










9,419, 776 


9,485,469 


9,534,279 


9, 544, 580 


9, 554, 571 


9, 568, 143 


286, 121 


j ivi ivi o Hj j \ 












9,422, 181 


9,501,327 


9, 553, 584 


9, 566, 004 


9,574,613 


9, 588, 975 


306,953 


(f( CL )) 


7 








8,445,057 


8, 570, 389 


8, 630, 159 


8, 674, 568 


8, 683, 940 


8, 693, 030 


8, 705, 378 


449, 167 


1 P\ Ivi IVI & Hi J \ 










8,448,582 


8, 576, 155 


8, 648, 195 


8, 695, 760 


8, 707, 065 


8,714,901 


8, 727, 973 


471, 761 




8 






8, 243, 496 


8,432,051 


8, 557, 190 


8,616,868 


8,661,208 


8, 670, 566 


8, 679, 642 


8,691,971 


1, 043, 242 










8, 229, 268 


8,421,009 


8,548, 167 


8,619,971 


8,667, 381 


8, 678, 649 


8, 686, 460 


8, 699, 489 


1, 050, 760 


(f( CL >) 


9 




8, 470, 989 


9, 129, 696 


9,338,521 


9,477, 113 


9, 543, 206 


9,592,313 


9, 602, 676 


9,612,728 


9, 626, 383 


3,950,814 


|^(MMSE)j 






8, 477, 596 


9, 121, 045 


9, 333, 566 


9, 474, 503 


9, 554, 088 


9, 606, 636 


9, 619, 125 


9, 627, 782 


9, 642, 223 


3, 966, 655 




1 i 


1.4925 


1.0778 


1.0229 


1.0148 


1.0070 


1.0051 


1.0011 


1.0010 


1.0014 




6,047, 061 






135.253 


33.803 


15.760 


19.847 


9.336 


2.001 


0.823 


0.219 


0.059 








-XMMSE) 
*j 


1.4937 


1.0759 


1.0233 


1.0151 


1.0084 


1.0055 


1.0013 


1.0009 


1.0015 




6, 139, 834 




(MMSE) 


132.917 


34.566 


14.742 


21.972 


8.547 


2.736 


0.789 


0.159 


0.061 







Table 2: Predicted parameter estimates, f , <r, cumulative chain ladder claims, Cjj, and estimated chain ladder reserves 
under the classical (CL) and Bayesian (MMSE) models. 



